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Abstract 

This paper proposes a dynamic congestion pricing model that takes into account mo- 
bile source emissions. We consider a tollable vehicular network where the users selfishly 
minimize their own travel costs, including travel time, early/late arrival penalties and 
tolls. On top of that, we assume that part of the network can be tolled by a central au- 
thority, whose objective is to minimize both total travel costs of users and total emission 
on a network-wide level. The model is formulated as a mathematical programming with 
equilibrium constraints (MPEC) problem and then reformulated as mathematical program- 
ming with complementarity constraints (MPCC). The MPCC is solved using a quadratic 
penalty-based gradient projection algorithm. 

1 Introduction 

According to the US Environmental Protection Agency (2006), in 2003, the transportation 
sector contributed to 27 percent of total U.S. greenhouse gas (GHG) emissions. This number 
is expected to grow rapidly with an estimated increase of transportation energy use by 48 
percent by 2015. Future transportation service network designs ought to take into account 
environmental issues. 

In this paper, we propose a dynamic second-best congestion toll problem with embedded 
emission model for the management and control of tollable vehicular networks. We assume 
that users of a given network are selfishly minimizing their own disutility, which consists of 
travel delay, early/late arrival penalties as well as the price of tolls. On top of that, there 
exists a central authority that undertakes the role of the Stackelberg leader, whose disutility 
includes two different aspects: the network efficiency and the automobile-induced emission. 

The upper-level decision variable for the central authority (Stackelberg leader) is a dy- 
namic congestion toll assigned to certain links of the network; while the lower-level decision 
variables for the travel agents (Stackelberg follower) include route and departure time choices. 
The proposed congestion pricing problem with embedded emission models is formulated as a 
mathematical programming with equilibrium constraints (MPEC) problem, with multiple ob- 
jectives including the mitigation of both congestion and traffic emission on a network-wide 
level. 
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To solve the multiobjective MPEC problem, we start by reformulating the differential 
variational inequality (DVI) - which is an equivalent form of DUE - into complementarity 
constraints. Then, we employ a weighted sum scalarization method to handle the multiple 
objectives. With these two steps, the multiobjective MPEC problem is transformed into a 
single-objective mathematical program with complementarity constraints (MPCC). To avoid 
the loss of constraint qualification, we relax the program by applying a quadratic penalty- 
based method. Then, the relaxed problem is solved with a gradient projection method in 
Friesz (2010). 

1.1 Congestion Toll Pricing 

The idea of employing toll pricing to mitigate both emission and congestion derives from 
congestion pricing strategy in Pigou (1920). In the literature, toll pricing problems can be 
classified into two categories: 1) First-best toll pricing, which tolls every arc in the network; 
and 2) Second-best toll pricing, which assumes that only a subset of the arcs is tolled for 
political or other reasons. Examples belonging to the first category were marginal social cost 
pricing strategy in Arnott and Small (1994), and other models and methodologies including 
Hearn and Ramana (1998), Dial (1999, 2000). As for the second-best tolling strategy, Law- 
phongpanich and Hearn (2004) proposed an MPEC approach to compute the optimal toll 
prices. However, the literature above focused on toll pricing in static models (see a compre- 
hensive review for static road pricing in Yang and Huang, 2005), which considered only the 
route choice decisions by the traffic agents, as noted in Yao et al. (2012). A more comprehen- 
sive and realistic model ought to take into account the time-dependent nature of traffic flow 
and dynamic toll pricing. Dynamic congestion pricing in bottleneck models were investigated 
in Arnott et al. (1990), Arnott and Kraus (1998) and Braid (1996), De Palma and Lindsey 
(2000). Based on link-delay model (Friesz et al. 1993, 2011), Friesz et al. (2007) proposed an 
MPEC model and a solution approach to determine the optimal second-best tolling strategy. 
Yao et al. (2012) further studied a dynamic congestion pricing problem with demand uncer- 
tainty. A more comprehensive review on the existing dynamic congestion pricing models and 
solution approaches was presented in Yao et al. (2012). 

The above literature provided evidences for the effectiveness of tolling strategies to stimu- 
late and optimize a traffic network. In this article, we seek to further exploit such a strategy to 
minimize both traffic congestion and automobile-induced emission simultaneously. To do this, 
we will provide an multiobjective MPEC model of a second best tolling strategy to determine 
the optimal toll price. Such an MPEC model incorporates the LWR-Lax model proposed in 
Friesz et al. (2012) to capture time-varying traffic behaviors. To the best of our knowledge, 
this is the first attempt to employ tolling to mitigate overall emission in the traffic network, 
the first dynamic toll pricing model that builds on LWR-Lax model in Friesz et al. (2012), 
and the first successful integration of LWR-Lax model with mobile source emission model. 

1.2 Dynamic User Equilibrium Model 

In this article, we employ the simultaneous route and departure time-choice DUE model pro- 
posed originally in Friesz et al. (1993). For the user equilibrium, unit travel cost, including 
early and late arrival penalties, is identical for those route and departure time choices se- 
lected by travelers between a given origin-destination pair. Such problem was articulated and 
formulated as a variational inequality (VI) in Firesz et al. (1993). The DUE model typi- 
cally consists of two major components: the VI formulation of equilibrium constraints and 
the network performance model known as the dynamic network loading (DNL) sub-problem. 
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The DNL aims at describing and predicting time evolution of system states by introducing 
dynamics to the traffic flows throughout the network. It determines the delay operator which 
is a key component of the DUE model. The DUE model is intertwined with the determination 
of path delay and thus the DNL. 

Friesz et al. (2001) studied the DUE as a single-level problem by formulating it as a 
controlled variational inequality. The DNL with link delay model (Friesz et al. 1993) and the 
route/departure time choice were solved simultaneously. In addition, the necessary conditions 
for optimal control problems with state-dependent time lags were derived. In Friesz and 
Mookherjee (2006), the theory of optimal control and the theory of infinite dimensional Vis 
were combined to create an implicit fixed point algorithm for calculating DUE. Friesz et al. 
(2011) extended the DUE problem from within-day time scale to a day-to-day time scale and 
formulated the dual time scale DUE problems. In Friesz et al. (2012), the Lighthill-Whitham- 
Richards model (Lighthill and Whitham 1955, Richards 1956) was employed in the DNL sub- 
problem, which was then formulated as a system of differential algebraic equations (DAEs) 
using the Lax-Hopf formula (Lax 1957, Evans 2010). The DUE problem was solved efficiently 
on reasonable size networks (Sioux Falls) with a realistic number of origin-destination pairs 
and paths. 

1.3 Mobile Source Emission Models 

Modeling approaches of mobile source emission can be classified into three categories: mi- 
croscopic, macroscopic and mesoscopic approaches. The microscopic emission models are 
relatively accurate: they characterize the emission rate on the level of a single vehicle, based 
on the physical attributes of the vehicle, driving behavior of the driver, as well as the sur- 
rounding environment. It is assumed that the rate of emission e(t) of a car can be expressed 
as a function of instantaneous velocity v (t) and acceleration a(t), 



such models can be easily calibrated and validated in the lab. There are several models based 
on microscopic emission (Barth et al 1996, Panis et al. 2006, Rakha et al. 2004). The 
drawback of the microscopic modeling approach is the lack of measurement associated with 
each individual car's dynamic. On the other hand, it is relatively easy to measure the traffic 
dynamics on a macroscopic level. The macroscopic emission models (Ekstrom et al. 2004) 
express the average emission rate e{f) on a road segment as a function of the average density 
p and average velocity v{t) in that same segment 



The drawback of the macroscopic modeling approaches for emission lies in the fact that the 
model is difficult to calibrate and validate, due to the lack of measurement of emission rate 
on a road. The third type of emission models - the mesoscopic emission model, approximates 
individual vehicles' dynamics using macroscopic flow models and measurements. Then the 
total emission rate is the aggregation of emission rates on a microscopic level, given by (jl.ip . 
The mesoscopic models (Csikos et al. 2011, Csikos and Varga 2011, Zegeye et al. 2010) take 
the modeling advantages of both macroscopic traffic flow models and microscopic emission 
models, thus overcome the drawbacks of the previous two approaches. 

In this paper, we will employ the mesoscopic point of view for emission estimation and 
prediction on a vehicular network. The emission estimation will be embedded in the Dynamic 
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Network Loading sub-problem, which describes and predicts macroscopic traffic variables given 
established path flows. The DNL provides a basis for the comparison of various microscopic 
emission functions, among which we distinguish between the two-argument functions e(t) = 
fi(v(t), a(t)) and the single-argument functions e(t) = /3(v(i)). 

The two- argument functions, such as the one proposed in the Modal emission model (Barth 
et al. 1996), use a physical approach that associates the power demand function to the var- 
ious driving conditions of the car such as low/high speed motion, acceleration, stop-and-go 
etc. The emission and fuel consumption functions are then expressed as affine functions of the 
power demand. Such model is very accurate, and can account for vehicles of different types. 
However, it is relatively difficult to integrate the modal model into macroscopic traffic flow 
models. In particular, the higher order traffic quantities such as acceleration cannot be suf- 
ficiently captured by first order models such as the Lighthill-Whitham-Richards conservation 
law (Lighthill and Whitham 1955, Richards 1956), see Section [3.11 for more discussion. 

The one-argument emission functions on the other hand, usually relies on the average 
speed. In Rose et al. (1965), it was shown when vehicle's travel speed is under 80km/hour, 
the relation between speed v (km/hour) and HC/CO emissions e x (pound/km) can be ap- 
proximated by (for now and sequel, e x denote the emission per unit distance). 

e x = hv- b2 (1.3) 

where b±, hi are parameters depending on vehicle type and surrounding environment. Kent 
and Mudford (1979) collected driving pattern data in Sydney and found that NO y emission 
e x can be modeled by 

~e x = h + ^ (1.4) 

v 

According to the Emission Factor Model 2000 (CARB 2000) by California Air Resources 
Board, constantly updated since 1988, the hot running emissions per unit distance 

e x = BERx exp|6i(t> - 17.03) + 6 2 (w - 17.03) 2 } (1.5) 

where BER stands for basic emission rates, which are constants associated with CO, NO y , HC. 
The unit of velocity is mile/hour, the unit of e x is gram/mile. 

1.4 Numerical Techniques for MPEC 

MPECs, with its bi-level nature, often create computational difficulty. A common approach 
to compute an MPEC is to formulate the bi-level program into a mathematical program with 
complementarity constraints (MPCC) (see examples in Ban et al., 2006, and Friesz, 2010). 
However, as noted in Rodrigues and Monteiro (2006) and in Ban et al. (2006), the comple- 
mentarity constraints might lose certain constraint qualifications. To resolve this issue, some 
regularization techniques were proposed. Ralph and Wright (2003) studied a relaxation ap- 
proach, which was applied by Ban et al. (2006) to solve a continuous network design problem. 
Anitescu (2000) proposed a Zi-penalty approach and studied its impact on the convergence 
of an interior point algorithm, while Monteiro and Meira (2011) tested a quadratic penalty 
function. According to their numerical results, quadratic penalty was a promising approach to 
handle complementarity constraints. However, discussions above all focused on MPCC in the 
context of finite dimensional programs. As for a continuous time dynamic MPEC, numerical 
techniques were scarcely visited, except for a single- level reformuation of Friesz et al. (2007), 
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a metaheuristic approach in Yao et al. (2012) and a simultaneous discretization-based method 
in Raghunathan (2004). 

In this paper, we will utilize the quadratic penalty method to solve our dynamic MPEC 
model. In specific, we will relax the model by dropping the complementarity constraints and 
attaching to the objective function a quadratic penalty function for the dropped constraints. 
Our numerical tests generate positive results. 

1.5 Organization 

The rest of this article is organized as follows. Section [2] briefly reviews the dynamic user 
equilibrium and its reformulation as variational inequality (Friesz et al. 1993) and differential 
variational inequality (Friesz and Mookherjee 2006, Friesz et al. 2010). The DNL sub-problem 
is based on a hydrodynamic network model proposed in Friesz et al. (2012), which is formu- 
lated as a system of differential algebraic equations (DAEs). In Section [3j two emission models 
are discussed in detail and embedded in the DNL sub-problem. In Section 2] and Section [SJ we 
will discuss the multi-objective MPECs, MPCCs and the solution method which is a quadratic 
penalty-based gradient projection algorithm. In SectionEl the solution of MPEC problem with 
multi-objectives on a small network is presented and illustrated. In particular, we will show 
the effectiveness of dynamic toll in reducing both congestion and emission. A Braess paradox 
is also observed. 

2 Dynamic User Equilibrium 

In this section, we will briefly review the DUE problem which serves as the lower-level com- 
ponent of our MPEC formulation. The proposed DUE model was formulated as a variational 
inequality (Friesz et al. 1993) and a differential variational inequality (Friesz and Mookherjee 
2006), then solved via a fixed-point algorithm in Hilbert space (Friesz et al. 2011). 

2.1 The DUE formulation 

Consider a planning horizon [to, tf] C 5?+. The most crucial ingredient of a dynamic user 
equilibrium model is the path delay operator, which provides the delay on any path p per unit 
of flow departing from the origin of that path; it is denoted by 

D p (t, h) for all p G V 

where V is the set of all paths employed by travelers, t denotes departure time, and h is a 
vector of departure rates. The path delay operators are obtained from the DNL procedure 
by summing up arc delays along the path. From these we construct effective unit path delay 
operators 

^ p (t,h) = D p {t,h) + F[t + D p (t,h) -T A ] forallpeP (2.6) 

where Ta is the desired arrival time. Introducing the fixed trip matrix {Qij : (i, j) S W) , 
where each Qij E !R + is the fixed travel demand, expressed as a volume, between origin- 
destination pair (i,j) € W and W is the set of all origin-destination pairs. Additionally, we 
will define the set V%j to be the subset of paths that connect origin-destination pair (i,j) € W. 

Denote the set of path flows h = {h p : p £ V}, we stipulate that path flows are square 
integrable: 

he[C 2 + [t , t f ]f l 
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We define the set of feasible flows by 

A = < h > : f f hp(t)dt = Qij for all G W I C [L\ [t ,t f ]) lVl (2.7) 

Let us also define the essential infimum of effective travel delays 

Vij = essinf [* p (t, h) : p £ Vij] for all (i, j) € W 

The following definition of dynamic user equilibrium was first articulated by Friesz et al. 
(1993). 

Definition 2.1. (Dynamic user equilibrium). .A vector of departure rates (path flows) 
h* G Ao is a dynamic user equilibrium if 

h* p (t) > 0,p G Vij => tf p [t, fc* (t)] = 

W^e denote this equilibrium by DUE , Aq, [to,tf]). 

Using measure theoretic arguments, Friesz et al. (1993) established that a dynamic user 
equilibrium is equivalent to the following variational inequality under suitable regularity con- 
ditions: 

find h* G A such that "j 

£ f tf %(t,h*)(h p -h* p )dt>0 \vi(%Ao,[to,tf]) (2.8) 
pevJto 

for all h G Aq J 

It has been noted in Friesz et al. (2011) that (|2.8p is equivalent to a differential variational 
inequality. This is most easily seen by noting that the flow conservation constraints may be 
re-stated as 

p £V *i } for all (i,j) G W 

2/y(*o) = 

Vij (if) = Qij ) 

which is recognized as a two point boundary value problem. As a consequence (|2.8p may be 
expressed as the following differential variational inequality (DVI): 



find h* 6 A such that 

ft 



/ / %{t, h*){h p - h* p )dt > I DVI(V, A, [t , t f ]) 
P ev Jt ° 



'to 

for all /i G A 



(2.9) 



where 



A = | /t > : ^ = K (*) > Wj(*o) = °> i/tf (*/) = for a11 (*'. J) e W [ ( 2 - 10 ) 

Finally, we are in a position to state a result that permits the solution of the DVI (|2.9p to 
be obtained by solving a fixed point problem: 
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Theorem 2.2. (Fixed point re-statement). Assume that *f> p (-,h) : [t OJ tt] — > K + is 
measurable for all p G V, h G A. Then the fixed point problem 

h = P A [h-a^(t,h)}, (2.11) 

is equivalent to DVI(^,A,A) where P\ [•] is the minimum norm projection onto A and a G 

Proof. See Friesz (2010). □ 
2.2 The DNL subproblem 

For the DNL sub-problem of the DUE model, we will employ the LWR-Lax model proposed 
in Friesz et al. (2012), which is a simplified version of the LWR model on networks. It is based 
on the assumption that any queues induced by congestion does not have physical size, thus no 
spill back occurs in the network. The network dynamics including route and departure time 
choices are summarized by a system of differential algebraic equations (DAEs). The derivation 
of the DAE system for DNL involves the variational method for Hamilton-Jacobi equations 
known as Lax-Hopf formula |Evans (2010) Lax (1957)] . For the brevity of our discussion, we 



will directly present the DAE system describing and predicting the flow propagation as well as 
path delay. More detailed derivation of this DAE system can be found in Friesz et al. (2012). 

Given a vehicular network (A, V), where A denotes the set of arcs, and V denotes the set 
of nodes. We define for each arc e G A, the free flow speed Vq and the jam density Pj am . 
Assume that the arc dynamic is governed by the following conservation law 

d tP e (t, x) + d x f{p e (t, x)) = (2.12) 

where p e (t, x) is the vehicle density on the link e. The fundamental diagram / e (-) is given as 
follows (Greenshields) 

f(p) = v e p(l--f-) (2.13) 



P 



jam 



We start by introducing some notations 



A : the set of all arcs 

V : the set of all nodes 

W : the set of all origin-destination pairs 

V : the set of utilized paths 

Vij : the set of utilized paths that connects origin-destination pair G W 

p = {ei, e2, • • • , e m (p)} € T 7 , e% G A : path represented by the set of arcs it uses 

h p (t) : departure rate at origin, associated with path p 

Qp(t) : entering vehicle count at arc e associated with path p 

qp(t) : enter rate at arc e associated with path p 

Wp(t) : exiting vehicle count at arc e associated with path p 

Wp(t) : exit rate at arc e associated with path p 

L e : length of arc e G A 
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ip e : the Legendre transformation of the following function 
<j)(u) = inf {p € [0, p e jam ] : f e (p) = u} 
where / e (-) is the fundamental diagram associated with e£i 
D(t;Q e ) : delay operator on arc e, associated with boundary profile Q e 

Q e {t) = w^(t) = y. w p^ 

By convention, we write qp^it) = h p (t), Wp°(t) = h p (t). T he following DAE system 
(|2.14p - (|2.19p for network loading was given in [Friesz et al. (2012)] . 

Q £ (t) = 52Q° p {t), q%t) = w e (t) = J2 W P® (2-14) 

egp e£p e£p 

j t Ql{t) = q;(t), jW e (t) = w e (t) for all peP (2.15) 

<#(t) = w?- 1 ®; ie[l,m(p)], peV (2.16) 

W e (t) = minjQ^^ + L 6 ^^^^)}; for all e£A (2.17) 

Q e (t) = W e (t + D(t; Q e ))- (2.18) 

t + D(t;Q«)) = ^-w^{t + D{t;Q e % i E [1, m(p)], pG? (2.19) 



We note that (|2.14p is definitional, i.e. the traffic on an arc is disaggregated by different 
route choices. (|2.16p represents the fundamental recursion, which allows the algorithm to 
carry forward to the next arc in the path. (|2.17p is the Lax-Hopf formula (Bressan and 
Han, 2011a,b). (|2.18p is often referred to as the flow propagation constraint, from which the 
travel time function D(-;Q e ) can be solved. (|2.19p describes the diverging model satisfying 
first-in- first-out . 

One shortcoming of the above DNL procedure is the lack of consideration for spillback, 
which not only aggravate congestion, but also cause a lot more stop-and-go waves (Colombo 
and Groli 2003) that will be another source of traffic emission. Thus as our next effort to 
model emission on a network level, we will consider the Lighthill-Whitham-Richards model 
on a network that captures spill over. 



3 The DNL subproblem integrated with emission models 

In this section, we will discuss two mesoscopic approaches for modeling traffic emission on a 
road network. The emission model will be considered in connection with the DNL sub-problem, 
which is consistent with established path flows as well as the flow propagation constraints. 
As a result, the output of the DNL sub-problem will include 1) the effective delay associated 
with each pair of departure time and route choices, and 2) the emission associated with each 
pair of departure time and route choices, as well as the total emission of the network. 



3.1 Emission as a functional of velocity and acceleration 

Consider a road network [A, V). For each arc a £ A, let us denote by p a (t, x), v a (t, x) the 
density and velocity of vehicles at time t and location x. The classical Lighthill-Whitham- 
Richards (LWR) model (Lighthill and Whitham 1955, Richards 1956) then describes the 
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temporal-spatial evolution of p a (t, x) by 

—p a (t, x) + —(p a (t, x)v( Pa (t, x))J = (3.20) 

where the velocity is expressed as an explicit function of density. The map p i— > p-v(p) is the 
fundamental diagram. 

Following the modal emission model on a microscopic level (Barth et al. 1996), the rate of 
emission e(t) of a moving vehicle can be modeled as a function of instantaneous velocity v(t) 
and acceleration a(t): 

e(t) = £(v{t), o(t)) (3.21) 

see Barth et al. (1996) for a detailed description of the function £. The total emission 
associated with a particular vehicle is expressed as a functional of velocity and acceleration 

ftf 

total emission = / £(v(t), a(t)J dt 

Consider an arc a £ A represented by a spatial interval [0, L a ] and a solution p a (t, x), v a (t, x), 
(t, x) E [Lq, tf] x [0, L a ] to the LWR conservation law (]3.20p . Then the total emission on this 
arc is given by 

p a (t, x) ■ e(t, x)dxdt (3.22) 



t Jo 



t JO 
tf i-L 



p a (t, x) ( £ ( v a (t, x), jjiv a (t, x) ) ) dxdt (3.23) 



= it l Pai 1 ^ x ) {^a, ^Va + Va- ^V^j dxdt (3.24) 

where = + v a ■ is the material derivative in Eulerian coordinates corresponding to 
the acceleration of the car in Lagrangian ones. 

Notice that the formulation (|3.24p can only be considered in the distributional sense, since 
the solutions p a , v a to the scalar conservation law are not differentiable. A detailed study of 
the continuous-time and discrete-time formulations of macroscopic emission models based on 
(|3.24p is left as future research. 



3.2 Emission as a functional of velocity 

In this section, we will discuss the travel speed emission models, and integrate it with the DNL 
sub-problem in Section [2.21 The model is based on the average travel speed for an arbitrary 
period of time. Such models ignore certain granularity of the flow dynamics and are calibrated 
and validated by empirical data, see Rose et al. (1965), Kent and Mudfor (1979) and CARB 
(2000) for more discussion on velocity-based emission models. We assume that the average 
emission rate e(t) is a function of average velocity v(t): 

e(t) = £{v(t)) (3.25) 

The speed v(t) can be averaged over a period that, say, spans a link traversal time. Given any 
feasible path flows ft £ A, we can solve the DNL solution in a way proposed in Friesz et al. 
(2012). Let D p (t, h) denote the total traversal time for path p given that departure from the 
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origin occurs at time t. In addition, we let Ta^(i) be the time of exit from arc given that 
departure from the origin occurs at time t and path p is followed, where p = {a\, . . . , a m u>)\- 
Then the average speed on the link ctj when the departure time from the origin occurs at t, 
denoted by v a .(t, h) is given by 

v ai (t,h) = P tts Lq ' p — 7-7, te[t ,t f ], per 

where L ai is the length of arc a% € p. Then using equality ()3.25j) . the contribution to emission 
of users departing at time t along path p is given by 

E r(t , H) = £ (rl (i) - <_,(*)) ■ g ( <[t) \ i[t) ) (3-26) 

Notice that we express the emission in a way similar to the effective delay operator ^ p (t, h), 
such formulation will facilitate the derivation of gradient of the objective function in Section 

El 

The total emission of the network associated with the path flow vector h is easily calculated 

as 



total emission = j h p (t) ■ E p (t, h) dt (3.27) 



Remark 3.1. In Rose et al. (1965), Kent and Mudford (1979) and CARB (2000), it was the 
emission per unit distance e x that was measured against travel speed, instead of emission per 
unit time et, see il.3\) . jl-4\) and fZHP- However, the emission per unit time as in \3.25\) can 
be easily calculated as 

— d dx d 

£(v(t)) = e(t) = Qf{t,x) = -qIq^<^ x ) = v{t)-e x (3.28) 
where e x is given by any one of hl.3\) . ^l-4\ ) or 111.5]) . 



4 Multiobjective Toll Pricing 

Most of the current MPEC-based Dynamic Traffic Assignment problems deal with a single 
objective. Lawphongpanich and Hearn (2004) studied the efficient tolls in a static network. 
Friesz et al. (2007) extended their work to consider dynamic congestion tolls and Yao et al. 
(2012) further studied a dynamic congestion pricing problem with demand uncertainty. In 
these studies, it was assumed that the Stackelberg leader seeks to minimize just the network 
congestion. However, the problem of congestion pricing with an emission-related objective is 
more subtle. The difficulty arises from the paradox that the most environment-friendly driving 
condition turns out to be very inefficient in terms of travel time (CARB 2000). Therefore, one 
major challenge we must overcome in solving the proposed model is to resolve the potential 
conflict between transportation efficiency and emission. Therefore, we proceed to formulate 
our dynamic congestion pricing problem as a multiobjective program: 



minU 

y 



to 



* P (t, h*)h* p {t)dt, f f Y E p(t, h*)h*(t)dt 



(4.29) 
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subject to 

rt 



£ [ \%(t, h*) + 5 a , p y a )(h* p -h p )<0 (4.30) 
he A (4.31) 



A := < 



h>0:^= J2 h P^> f«(*o) = °. W (*/) = <5y for a11 G W f ( 4 - 32 ) 

< y« < Y UB for all a e ^ (4.33) 



where y := (3^ a : a £ ^4), the first term appearing on the right hand side of (|4.29p is the total 
effective delay, the second term is the total emission as given in (|3.27j) . and where 



1 if path p traverse arc a 
otherwise 



and Yjjb £ ^+ denotes the prescribed upper bound to the toll. Constraint (|4.30j) is a straight- 
forward extension of the VI formulation of DUE, which takes into account the toll prices y a . 
One crucial term in the formulation is the effective delay operator \&(-, •) : A — > (C[to, tf]y, '. 
The evaluation of such operator is referred to as the DNL procedure and is realized by solv- 
ing the DAE system in Section 12.21 Therefore the above program is further constrained by 
(|2.14p - (|2.19p . and (|3.26p - (13.27p . Notice that in the program above, U is a vector of objective 
functions. 

Definition 4.1. (Pareto Optimal) For a multiobjective optimization problem of the form: 

minF(x) = [F 1 (x),F 2 (x),...,F k (x)f 

subject to 

a feasible solution x* E X , where X denotes the feasible region, is Pareto optimal if and only 
if there does not exist another solution, x G X , such that F{x) < F{x*), and Fi(x) < Fi(x*) 
for at least one function. 

A Pareto optimum requires that no other feasible solutions that improve at least one 
objective without deteriorating another. Seeking to attain a Pareto optimum, we employ a 
common approach called weighted sum method. Successful examples of such an approach can 
be found in Zadeh (1963) and Murata et al. (1996). Some new insights into the weighted sum 
scalarization can be found in Marler and Arora (2010). 

5 Solution Methodology 

The DVI can be reformulated as complementarity constraints as follows: 

(%(t, h*) + 5 aiP y a - fUj) ±h; for all p G Vij, ij G W (5.34) 

%(t, h*) + 5 a , p y a - fXij > for all p G Vij, ij G W (5.35) 
h* p > for all p G V (5.36) 
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where h* € A. With complementarity constraints substituting the DVI in the MPEC model, 
we are able to obtain a single level mathematical program consisting of the objective function 
(|4,29p and (|2.14|) - (|2.19|) . and (|3.26|) - (|3.27|) . Since the complementarity constraints may not 
satisfy MFCQ (Rodrigues and Monteiro, 2006, and Izmailov and Solodov, 2004), we propose 
to apply a penalty method to handle these constraints. In Monteiro and Meira (2011) the 
quadratic penalty approach, or sometimes called sequential penalty technique, was tested 
numerically with positive results in solving an MPCC. Following the quadratic penalty method, 
we penalize the complementarity constraints and obtain an augmented objective function as: 



where 



U := \XJ x {h\y,^M), U 2 (h*,y,fi,M)] (5.37) 

U^h\y,^ t M) := £ Y / y p (t,h*)h* p dt + Q(h*,y,ti,M) (5.38) 
ijeWpeVij *° 

U 2 (h*,y,^M) := Y I S E p {t,h*)h* p dt+Q(h*,y,n,M) (5.39) 



Q(h*,y,n,M) :=M Y Y P ' [OMW + S ajP y a - 2 dt 

+ m Yj y p ^ - w - Sa <p y *> dt ( 5 - 4 °) 

where [i := (/Xy : ij G W), and M is a properly large number. In order to compute the 
above multiobjective problem, we use a simple but commonly-used weighted sum scalarization 
method: 



s u (h,y, l M,M) = aU 1 (h*,y,(ji,M) + pu 2 (h*,y,n,M) (5.41) 

where a, f3 £ are weights for two objectives respectively. For normalization, we further 
require that a + f3 = 1. Then, the problem becomes a single-level single-objective optimal 
control problem. Such a problem is computed with the gradient projection method in Friesz 
(2010). 



6 Numerical study 

In this section, we will numerically illustrate the solutions to the MPEC problems as well as 
the effectiveness of the dynamic toll in mitigating both congestion and emission. Consider the 
network shown in Figure [lj which consists of six arcs and five nodes. There are two origin- 
destination pairs (1, 3), (2, 3), among which four and two paths are utilized, respectively. 

Pi = {3, 6}, p 2 = {1, 2, 6}, p 3 = {1, 2, 4, 5}, p A = {3, 4, 5}, p 5 = {6}, p 6 = {4, 5} 

Assume that arc 1 is tolled, the upper-level decision variable of the MPEC problem is the 
dynamic toll price imposed on arc 1. The lower level is a DUE problem where drivers choose 
their departure time and route in order to minimize the total travel cost, including a toll price. 
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Figure 1: The 6 arc 5 node network 



6.1 Numerical setting 

We fix a morning commute horizon spanning five hours from 6:00 am to 11:00 am. The arc 
parameters are shown in Table 1. 



Arc 


Jam density 


Free flow speed 


Length 




(vehicle/mile) 


(mile/hour) 


(mile) 


1 


400 


35 


10 


2 


400 


35 


10 


3 


400 


35 


10 


4 


400 


35 


20 


5 


400 


35 


20 


6 


400 


35 


15 



Table 1: Arc parameters for the network 



We will employ the emission model discussed in Section 13.21 and Remark 13.11 

£(v(t)) = v{t)-e x 

where the hot running emission e x is given by (jl.5[) : 

e x = BER x exp{6i(u - 17.03) + 6 2 (> - 17.03) 2 } (6.42) 

where BER = 2.5, b 1 = -0.04, b 2 = 0.001 (Smit 2006). We consider two cases in our 
computation: 

I. The demand matrix is (Qi,3, ^2,3) = (820, 410), and The upper bound for toll price 
is Y UB = 10. 

II. The demand matrix is (Qi,3, Q2,3) = (1400, 700), and The upper bound for toll price 
is Y UB = 10. 

6.2 Numerical results 

The solution algorithm for MPCC is implemented in Matlab (2010a), which runs on the 
Intel Xeon 3160 Dual-Core 3.0 GHz processor provided by the Penn State High Performance 
Computing center. 
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6.2.1 Case I 



We display the first numerical result from Case I above in Figure El [21 [3] and [5j For comparison 
reasons, we first solve the DUE problem without any tolls. The departure profiles at the 
beginning of each path are shown in Figure [2j Then the optimal solution to the MPEC 
problem which is a dynamic toll on arc 1 is presented in Figure 01 Under such toll price, 
the new equilibrium solution is now displayed in Figure O Notice that two paths p2 , P3 use 
link 1 thus are affected by the toll. In the MPEC solution with toll, the path flows on p2, pz 
diminish to the point where path p% is not used and very little traffic uses path p2- Figure 
[5] shows the difference between the equilibrium path flows without toll, and the equilibrium 
path flows with optimal toll. It is clearly that most of the users on path p2, ps switch to path 
pi, P4, due to the presence of toll. 

We also compare the two objective functions in the case of DUE with toll and DUE without 
toll. The results are summarized in Table 2. It is clear that by adding the toll, we are able to 
reduce the total travel cost and total emission by 2.9% and 10.4% respectively. 





Figure 2: Case I. DUE solution without any Figure 3: Case I. DUE solution with optimal 
toll toll. 




Path 1 






Path 2 






Path 3 






Path 4 






Path 5 






— Path 6 




1. 


I 






U7 



1.5 2 2.5 

Time (hour) 



Figure 4: Case I. Optimal toll on arc 1. 



Figure 5: Case I. Differences of path flows be- 
tween DUE without toll and DUE with toll. 
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Total travel cost 


Total emission 


DUE without toll 


3.4744E+04 


3.1789E+06 


DUE with toll 


3.3723E+04 


2.8483E+06 



Table 2: Comparison between DUE solutions. Case I. 



6.2.2 Case II 

In the Case II, the demand for each O-D pair is increased. The same quantities of the solution 
as in Case I are shown in Figure El [8] and EE Unlike the first case, in Case II there is 
very little difference of the DUE solutions with and without tolls. We interpret such results 
with the following intuition: when the demand increases, the system becomes less sensitive 
to control parameters, making the system less controllable. This is also reflected from the 
comparison of objectives, as shown in Table 3. The reduction of total travel cost and total 
emission is only 0.04% and 0.45%. 



- Path 1 

- Path 2 

- Path 3 

- Path 4 

- Path 5 
■ Path 





Figure 6: Case II. DUE solution without any Figure 7: Case II. DUE solution with optimal 
toll. toll. 
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Figure 8: Case II. Optimal toll on arc 1. 



Figure 9: Case II. Differences of path flows 
between DUE without toll and DUE with toll. 
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Total travel cost 


Total emission 


DUE without toll 


7.5962E+04 


5.4119E+06 


DUE with toll 


7.5932E+04 


5.3878E+06 



Table 3: Comparison between DUE solutions. Case II. 



6.3 Different weights, the Braess's paradox 

The multi-objective program is solved using the weighted sum scalarization method. We use 
different weights for total travel cost and total emission and see how the solution is affected. 
The test is conducted for Case I and Case II, the results are summarized in Table 4 and Table 
5, respectively. We indicate by a the weight for the total travel time, by /3 the weight for the 
emission. 





Total travel cost 


Total emission 


a 


P 


Weight i, 


3.3723E+04 


2.8482E+06 


0.0988 


0.9011 


Weight ii, 


3.3723E+04 


2.8483E+06 


0.9434 


0.0566 



Table 4: Comparison of objectives for different choices of weights, Case I. a is the weight of total 
travel cost, /3 is the weight of total emission. 





Total travel cost 


Total emission 


a 


P 


Weight i, 


7.5932E+04 


5.3878E+06 


0.0138 


0.9862 


Weight ii, 


7.8360E+04 


5.2396E+06 





1 


Weight iii, 


7.5858E+04 


5.4175E+06 


1 






Table 5: Comparison of objectives for different choices of weights, Case II. a is the weight of total 
travel cost, (i is the weight of total emission. 

Another interesting phenomenon is the Braess-like paradox displayed in Case I, where the 
performance of the network, whether in terms of travel cost or emission, is enhanced in the 
more constrained system (the one with toll). 

6.4 Discussions 

In the numerical tests, we solved two cases of a simple network. The results provide some 
interesting phenomena and interpretations to our proposed model and the dynamic toll pricing 
problem in general. 

First, as noted above, the optimal toll creates a Braessdike paradox in Case I. A Braess' 
paradox (Braess 1969) states that the system overall performance can be reduced in some 
cases with added extra capacity when traffic agents selfishly decide their route. In our problem 
setting, the traffic agents all choose both their route and departure time that appear most 
favorable to themselves. By imposing the toll, affordable time windows on path p2 and path 
P3 become significantly smaller. In particular, we notice that in the presence of the toll, the 
path j>3 is completely abandoned by the users. With fewer affordable choices on path and/or 
departure time, the system's capacity becomes more constrained, yet two measures of the 
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system performance, the total traffic cost and the overall emission amount, are reduced, i.e. 
a higher efficiency is attained in terms of both transportation and environment. 

Second, tolls can act as effective stimuli in a transportation system. We observe from 
Table 2 (in Case I) and Table 3 (in Case II) that by properly choosing the toll prices, one can 
reduce both the total traffic cost and total emission. In Case I, the price of toll on the first 
arc of path P2 renders this path no longer an affordable choice, i.e. the equilibrium flow on 
this path vanishes, which, nonetheless, creates a Braess-like paradox as discussed above. 

Finally, by comparing Table 4 and Table 5, in response to changing weights for the two 
objectives in our weighted sum approach, Case I (in Table 4) shows very minor changes in 
objective values compared to Case II (in Table 5). The reason for such a difference in the 
sensitivity to weight perturbations is still of our research curiosity. Nonetheless, this points to 
a necessity of a careful determination of weights for the two objectives in our model in order 
to attain the optimization goal of the central control. 

7 Conclusion 

This paper proposes a congestion pricing problem that takes into account the environmental 
impact of traffic dynamics on a vehicular network. An MPEC problem is formulated which 
aims at mitigating both emission and congestion on a network level. For the vehicle-induced 
emissions, we employ a mesoscopic modeling approach which uses macroscopic traffic quanti- 
ties and microscopic emission functions. The emission model is then embedded in the Dynamic 
Network Loading sub-problem of DUE in the following way: for every given set of departure 
rates at the origin, we run the network loading procedure until all macroscopic traffic quanti- 
ties are available; then we provide information of vehicle dynamics to the microscopic emission 
function and perform emission estimation on a network-wide level. The MPEC problem is 
first reformulated as an MPCC. To avoid violation of constraint qualifications, we then ap- 
ply a quadratic penalty-based method to relax the original program. The relaxed model is 
solved with the gradient projection algorithm (Friesz 2010) with multiobjective handled via 
a weighted sum scalarization. The numerical examples have demonstrated the effectiveness 
of congestion toll in controlling and reducing both total travel cost and emission. We also 
observe an instance of the Braess-like paradox where a more constrained system results in 
higher transportation efficiency and less environmental impact. 

The emission model employs the average velocity-based emission function (CARB 2000). 
Such model is computationally convenient in terms of integration with the DNL procedure. 
However, it relies on the average speed of cars, therefore is less accurate compared to other 
models that depends on instantaneous speed and acceleration. As our future research, we will 
take into account spill over into our DNL model and capture the speed variation in a more 
accurate way. 

Kumar and Vladimirsky (2010) pointed out that a weighted sum scalarization approach 
could obtain only convex part of the Pareto front, which might lead to selecting suboptimal 
trajectories. To resolve such an issue, they provided an alternative 'marching' method. The 
application of their approach to our model is also of our future research interest. 
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